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We present explicit kinetic equations for quantum transport through a general molecular quantum- 
dot, accounting for all contributions up to 4th order perturbation theory in the tunneling Hamilto- 
nian and the complete molecular density matrix. Such a full treatment describes not only sequential, 
cotunneling and pair tunneling, but also contains terms contributing to renormalization of the molec- 
ular resonances as well as their broadening. Due to the latter all terms in the perturbation expansion 
are automatically well-defined for any set of system parameters, no divergences occur and no by- 
hand regularization is required. Additionally we show that, in contrast to 2nd order perturbation 
theory, in 4th order it is essential to account for quantum coherence between non- degenerate states, 
entering the theory through the non-diagonal elements of the density matrix. As a first application, 
we study a single-molecule transistor coupled to a localized vibrational mode (Anderson-Holstein 
model). We find that cotunneling-assisted sequential tunneling processes involving the vibration 
give rise to current peaks i.e. negative differential conductance in the Coulomb-blockade regime. 
Such peaks occur in the cross-over to strong electron-vibration coupling, where inelastic cotunneling 
competes with Franck-Condon suppressed sequential tunneling, and thereby indicate the strength 
of the electron-vibration coupling. The peaks depend sensitively on the coupling to a dissipative 
bath, thus providing also an experimental probe of the Q-factor of the vibrational motion. 



PACS numbers: 73.63.Kv , 85.65.+h , 63.22.-m 
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I. INTRODUCTION. 



Electron transport through single-molecule transistors 
(SMTs) has been intensively studied theoretically in re- 
cent years [J E S, 0, H, i, 0, EM ^M^L^ivenby 
ongoing experimental advances[13l. Il4 Hal 1 6l. 1 1 7l. Il8l. Illf . 
One of the most distinctive features of SMTs, compared 
to artificially nano-structured devices such as quantum 
dots, is the coupling between their quantized mechan- 
ical and electronic degrees of freedom The size 
and shape distortions of an SMT [l7| and its center-of- 
mass motion [lH result in sharp transport resonances 
whose amplitudes are governed by the quantum mechan- 
ical overlap of the corresponding mechanical wavefunc- 
tions. This Franck-Condon (FC) transport effect is of 
fundamental interest since it is induced by the change 
of molecular charge, therefore involving strong electron 
charging and non-equilibrium effects, in contrast to the 
usual FC-effect in optical spectroscopy where the charge 
remains unaltered. The discrete vibrational modes of 
a molecule are also important in assessing the atom- 
istic details of the transport junction [2fj| . Finally, the 
demonstrated control over the molecular energy levels of 
an SMT using a gate electrode provides interesting per- 
spectives for realizingquantized nano-electromechanical 
systems (NEMS) HUES @l ■ 

The basic FC transport picture [24| assumes single 
electrons to sequentially tunnel on and off the SMT. This 
is valid in the limit of weak tunnel coupling and for ap- 
plied gate- and bias- voltages such that fluctuations of the 



molecular charge are not suppressed by Coulomb interac- 
tion (Coulomb blockade) or quantum confinement effects. 
In this limit it is sufficient to describe transport in low- 
est non-vanishing order perturbation theory in the tun- 
neling and many interesting results have been reported. 
For instance, a well studied model in this context is the 
Anderson-Holstein model, consisting of a spin-degenerate 
level with a linear coupling (electron-vibration coupling, 
A) between the charge on the level and the coordinate of 
a vibrational mode. When the overlap integrals between 
low lying vibrational states in two adjacent charge states 
of the SMT vanish, a suppression of single-electron tun- 
neling (SET) occurs, called Franck-Condon blockade 0]. 
Here electron transport was found to take place- through 
self-similar avalanches, leading to bunching of electrons 
and enhanced shot noise. Extending the basic model with 
a charge-dependent vibrational frequency, additional res- 
onances occur Q, and interference of vibrational wave- 
functions was shown to lead to a suppression of the elec- 
tric current at finite bias [l[ due to a population inversion 
of the vibrational distribution. More complex models 
with (quasi-)degenerate electronic orbitals and multiple 
modes exhibit (pseudo-)Jahn- Teller physics. These may 
show rectifying behavior Q , dynamical symmetry break- 
and current suppression due to Berry phase ef- 
5|. Finally, distinctive transport signatures of the 
breakdown of the Born-Oppenheimer separation [|| and 
correlations of vibration- and spin-properties have been 
predicted, such as a vibration-induced spin-blockade [i"2| . 

Since the complicated transport processes in SMTs 
may result in a suppression of single-electron tunnel- 
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ing, it becomes more urgent to understand the effect 
of higher order tunnel processes. Even more so since 
experimentally SMTs typically exhibit a significant tun- 
nel coupling. The purpose of this paper is to set up a 
general method to properly describe all coherent tunnel 
processes in leading and next-to-leading order in the tun- 
neling Hamiltonian. The method applies to very general 
molecular quantum dot models, with many quantized ex- 
citations and few relevant selection rules for transport 
quantities. Additionally, all molecular interactions in- 
cluded in the model, for instance Coulomb charging and 
electron-vibration, are accounted for exactly by from the 
outset formulating the transport equations in the basis 
of many-body eigenstates of the molecular Hamiltonian. 
The non-linear transport is obtained using the explicitly 
calculated non-equilibrium density matrix. A variety of 
next-to-leading order effects have been discussed previ- 
ously. For instance, a well known signature of higher 
order tunneling processes is the appearance of inelas- 
tic cotunneling steps in the differential conductance, the 
position of which arc independent of the gate voltage, 
which have been observed in many experiments on semi- 
conductor and molecular quantum-dots [HI, 0, OH, [HI . 
Additional gate-voltage dependent transport resonances 
have been found inside the Coulomb blockade regime (25[ 
and discussed theoretically [26|, [2?| . These resonances are 
due to sequential tunneling events starting from states 
excited by inelastic cotunneling processes ("heating the 
molecule" ), called cotunneling- assisted sequential tunnel- 
ing (COSET). It was recently suggested [28[ that these 
resonances can be used to experimentally estimate the 
decay-rate of vibrational excitations due to a coupling 
to a dissipative environment. Indeed, by extending the 
Golden-Rule approach by a next-to-leading order expan- 
sion of the T-matrix in the tunneling, it was found Q 
that the COSET features are particularly pronounced 
in the Andcrson-Holstcin model in the limit of large 
electron-vibration coupling. Finally, effects of electron 
pair-tunneling were discussed [Io| for an effective Ander- 
son model with attractive electron-electron interaction in 
the Golden-Rule approach using a Schrieffer- Wolff trans- 
formation. 

The method set up here captures all the aforemen- 
tioned effects simultaneously. By providing a microscopic 
derivation we overcome some drawbacks of the meth- 
ods employed in the cited works, related to account- 
ing for broadening and renormalization of the molec- 
ular resonances, which were previously discussed in 
e.g. @, IH, HO, [3l|], see also Ref. [H. The main focus 
of the paper is therefore on the general aspects of the 
transport theory. As a central result we present explicit 
kinetic equations from which the full molecular density 
matrix and transport current can be calculated. We re- 
formulate the real-time transport theory [H, HH using 
Liouville super-operators to present a straight-forward 
derivation. The expressions for the transport rates are 
valid for a wide class of quantum-dot systems and, im- 
portantly, involve no assumptions on model-specific se- 



lection rules. We show that in such higher order calcula- 
tions it is crucial to include contributions from coherent 
superpositions of molecular states not protected by se- 
lection rules, even when the level spacing is much larger 
than the tunneling broadening. The Anderson-Holstein 
model with a large vibrational frequency compared to the 
tunneling coupling, fku T, presents a case where this is 
extremely important and we demonstrate our method for 
this model. This is in clear contrast to lowest order per- 
turbation theory where these so-called non-secular terms 
can be neglected. We are not aware of previous works 
pointing this out. 

To maintain readability, the paper is divided into three 
parts: a general part [H] application IIIII and technical 
Appendices. In Sect. [IT] we shortly describe the general 
model of a molecular quantum dot system and the ba- 
sic equations of the real-time transport theory. We then 
discuss the central results, the explicit transport equa- 
tions for the full density matrix and transport current. 
Detailed derivations and expressions are presented in a 
coherent way in the Appendices for the theoretically in- 
terested reader. In Sect. IIIII we study in detail the specific 
model of a molecular transistor coupled to a localized vi- 
brational mode, the non-equilibrium Anderson-Holstein 
model. We summarize and conclude in Sect. ITVl 

Throughout the rest of the paper we use natural units 
where h = &b = |e| = 1 where — |e| is the electron charge. 

II. MODEL AND TRANSPORT THEORY 

We consider a molecule as a complex quantum dot, 
connected to a number of macroscopic reservoirs labeled 
by r. The electrons in the reservoirs are considered to 
be non-interacting, but no assumptions are made con- 
cerning the type or strength of the interactions on the 
molecule, as long as we can diagonalize the isolated 
molecular many-body Hamiltonian. The entire system 
is described by the Hamiltonian H to t = H + Hr + Ht, 
where = J2 r and 

H = J2 E aW)(a\, (1) 

a 

H r = E / dw wc rCT _ w c r<T+w , (2) 

a'e{N- V ) 

Ht = E E E W dujT r %\a)(a'\c ra ^.(3) 

raNr 1 =±l a£N ^ 

The Hamiltonians are written from the outset in a form 
which deviates from that commonly used. This allows 
crucial simplifications of the derivations and explicit ex- 
pressions presented in the Appendices. In the molecular 
Hamiltonian, H, \a) denotes a general many-body eigen- 
state with energy E a . We assume that we can classify 
these states by the number of excess electrons, N . on 
the molecule. The electron number, together with other 
quantum numbers depending on the model at hand (e.g. 
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spin, magnetic and vibrational quantum numbers), are 
labeled by a. We will loosely denote by N a the electron 
number in state a. H r describes reservoir r and is written 
in terms of continuum field operators 



t+u = -p S(e rak - u)Crak, V = +, (4) 

k y pra 

r-o> = J2^= s ( er,jk _w ) c Lfc> v = -, (5) 



where c ra . k (c ra k) are the usual creation (annihilation) 
operators for electrons in reservoir r with spin-projection 
<7, state-index k and energy e r(7 fc. We will refer to rj as 
the electron-hole (e-h) index. p ra is the density of states. 
Inserting ([4][5]) into ([2]) one recovers the standard form 
of the reservoir Hamiltonian. For this one assumes that 
there is a unique correspondence between k and e r ak- 
For cases where this does not hold, one labels different 
branches of the dispersion relation by an additional in- 
dex. The tunnel Hamiltonian describes the tunnel- 
ing into or out of the molecule, involving a change of the 
molecular state from a 1 to a. The relevant matrix ele- 
ments are given by superpositions of single-particle tun- 
neling matrix elements ti ra and many-body amplitudes 
of the molecular wavefunction: 

T r a :'+ = V/^E^aMlla'), (6) 
l 

T^L = VP^T, *?ra (a\diM) = * • (7) 



Here I labels a single particle basis for the molecule with 
corresponding creation / annihilation operator d\ , di a . 
Note that the density of states is incorporated in T^ v , 
simplifying many^ expressions. The spectral densities 
r ab,aV _ 2nT™' + T^, thus represent the set of rele- 
vant energy scales for the tunneling. Both p ra and ti rr7 
are assumed to be energy independent. This is the most 
relevant physical limit and presents no principle limita- 
tion of the presented method (only numerical). Charge 
conservation implies the selection rule, N a — N a > = 1, 
which is contained in (a\d\ a \a') cx <5jv q ,jv This is the 
only selection rule assumed. A Fermion sign, ?y, appears 
in Eq. ^ since we always write the reservoir operator to 
the right of the projector. However, one can show that 
this exactly cancels in all expressions involving an aver- 
age over the reservoir degrees of freedom. It cancels with 
an extra Fermion sign appearing when disentangling the 
dot and reservoir operators, using that an equal number 
of creation and annihilation operators must occur to give 
a non-zero average, see Ref. [35| for a proof. We can 
therefore discard the sign 77 from the outset, and every- 
where treat dot and lead operators as commuting, greatly 
simplifying the calculation of signs. 



A. Kinetic Equation 

A microscopic molecular system coupled to macro- 
scopic reservoirs is completely described by its reduced 
density operator P(t), obtained by averaging the total 
density operator over the reservoir degrees of freedom, 
P{t) = TrRp(i). The reduced density operator evolves in 
time according to a quantum kinetic equation. The pres- 
ence of strong non-equilibrium effects (non-linear trans- 
port) and strong local interactions (Coulomb, electron- 
vibration, etc.) makes the calculation of the transport 
rates occurring in this equation a cumbersome task. 
Here our goal is to derive explicit expressions for the 
next-to-leading order transport rates in terms of the pa- 
rameters E a ,T^" of the Hamiltonians (Q][3]) and the 
statistical properties of the electrodes T (temperature) 
and p r (chemical potential) . The real-time transport 
theory, developed in [33|, Hj] and extended by several 
groups [U Hy, provides straightforward rules for 
the calculation of the transport rates using a diagram- 
matic representation on the Keldysh contour, avoiding 
any spurious regularization problems. This technique has 
been simplified further by using special Liouville super- 
operators and corresponding diagrams in the context of a 
non-equilibrium renormalization group approach [35l. l38l] . 
For clarity we discuss the general aspects here, whereas 
the important but cumbersome expressions are coher- 
ently derived and presented in the Appendices. The 
starting point is the time evolution of the density op- 
erator of the total system, molecule + reservoirs: 



P(t) 



— iLtatt 



p(0). 



(8) 



Here the Liouvillian super-operators in L tot — L + Lr + 
Lt, act on an arbitrary operator A by forming a commu- 
tator with the Hamiltonian, e.g. LA = [H,A\. We as- 
sume the system to be decoupled at the initial time t = 0, 
such that the density operator factorizes, p(0) = P(0) j or 
where pj\ = J7 r p r and p r describes reservoir r. Each 
reservoir is assumed to remain in internal equilibrium 
independently and is described by a grand-canonical en- 
semble at all times. When a bias voltage is applied to 
the system, causing the chemical potentials of different 
leads to differ, this puts an inhomogeneous "boundary 
condition" on the molecular density operator and drives 
it out of equilibrium. We now take the Laplace transform 
of Eq. ijSJ) and trace out the reservoirs 

/•oo 

P{z) = Tr / dte lzt e- lLtott P(0)p n (9) 
R Jo 

-P(0) PR (10) 



iTr 

R z — Lr — L — L\ 

P(0), 



z- L-iW(z) 



(11) 



where the last expression is obtained by expanding the 
denominator in (jTUJ) in powers of the tunneling Liouvil- 
lian, Lt, carrying out the trace over the reservoirs and 
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re-summing the series, see Appendix lAl for details. Here 
iW(z) is a (super-operator) self-energy and L + iW(z) 
describes the molecular density operator in the presence 
of the reservoirs. If Eq. (fTTj) is transformed back to the 
time-domain, W(t—t') appears as a kernel in the integro- 
diffcrcntial equation for P(t) 



P{t) = -iLPlt) 



dt'W(t-t')P(t'), 



(12) 



assuming that the Laplace transform of W(t — t') exists. 
We are exclusively interested in the stationary state at 
t — > oo (asymptotic solution) of the molecular density 
operator, i.e. the zero frequency limit z — > iO where 
the imaginary infinitesimal physically originates from the 
adiabatic switching on of the tunneling. Assuming that 
a unique stationary state exists and using lim t _ (00 F(i) = 
— i\inu_ fio zP(z), Eq. pip gives the standard form [30l . 
[331 - [35j of the stationary state equation: 



= (-iL + W)P. 



(13) 



Here, and in the rest of the paper, we use the notation 
W = \im z ^ioW(z) and P = limt_ >00 P(t). Supplemented 
with the probability normalization condition Trjyj-P = L 
where Trjyi is the trace over the molecular degrees of free- 
dom, this uniquely determines the stationary state. The 
normalizability derives from the general property of the 
kernel TrjyiVFA = for any operator A. Matrix elements 
of a super-operator, S, are defined according to 



S, 



a'b' 
a b 



(a\(S\a')(b'\)\b), 



(14) 



meaning that we first act with S on a projector |a')(6'|, 
generating a new operator, and subsequently take matrix 
elements of this. In the basis of the many-body eigen- 
states of the isolated molecule the molecular Liouvillian 



- a'b' 



[E a — Eb)5 aa i5bb< 



(15) 



Our main objective is to calculate the expectation value 
of the electron current flowing out of reservoir r into the 
molecule, I r (t) = Trl r p(t) where Tr is the trace over 
the full system and I r = —^N r = —i[H T ,N r ], with N r 
being the number operator for electrons in reservoir r. 
As is shown in Appendix [AJ this expectation value can 
be obtained from a kernel similar to W and the density 
operator. In the stationary state: I r = Tim{Wi t P}, 
where the current kernel, Wj T , contains the subset of 
tunneling processes described by W which contribute to 
the current through reservoir r. We can now write down 
the generalized, formally exact, master equations 




1 

Ir 



E 



— \ /a 



k=l 



a b 
b 



V P 

/ j 1 aai 



EEEK 



a a'b' k 



,(2*0 



a'b 



Pa'b' 



Pa'b', (16) 

(17) 
(18) 



Here we have expanded the kernels in even order terms 
2k in the tunneling Liouvillian accounting for coherent k- 
clcctron tunnel processes (odd orders vanish when trac- 
ing over the reservoirs since the tunneling Hamiltonian 
is linear in reservoir field operators). Eq. (p~6l Tl8]) com- 
pactly formulate the transport problem. The first central 
result of this paper is the explicit evaluation of the ker- 
nels and as given in Appendix (Eq. ([C3])) 
and [Dl (Eq . (|D4[) L accounting for coherent single- and 
two-electron tunneling processes. Their detailed form is 
not needed here, but an important property of these ex- 
pressions is that they are finite by construction for any 
system parameters and applied voltages at non-zero tem- 
perature: no divergences occur and no by-hand regular- 
ization is required at any stage of the calculation as is 
the case in the Golden- Rule T-matrix approach Q. Of 
course, the finite temperature must be chosen sufficiently 
large compared to the tunneling couplings (r <C T) to 
avoid the breakdown of perturbation theory. 

For the solution of the kinetic equation it is impor- 
tant to know whether the molecular density matrix is 
diagonal in certain quantum numbers due to a con- 
servation law. The only such law explicitly enforced 
here concerns the total charge in reservoirs + molecule: 
[i?tot; Ntot] = 0. As is shown in Appendix [Bl the matrix 
elements of W£ b b vanish unless the charge differences are 
equal: N a ' — Ny = N a — Nb- With the assumption that 
the density matrix is diagonal with respect to charge at 
t = 0, before the coupling to the reservoirs is switched 
on, it is guaranteed to remain so at all times. In a similar 
way any conserved quantity of the total system encodes 
selection rules in the tunneling matrix elements ensuring 
that the density matrix remains diagonal in the corre- 
sponding molecular quantum number. For example, for 
the Anderson-Holstein model studied in Scct. lIIII the con- 
servation of total spin-projection, S z , leads to a density 
matrix which is diagonal in the spin-projection of the 
molecule, s,. 



B. Solution of the Kinetic Equation 



The solution of equations ljl6M18[) with perturbativcly 
calculated kernels (up to a finite order) for the full molec- 
ular density matrix requires some care for models with 
excited states and tunnel matrix elements without strict 
selection rules. The present section is therefore devoted 
to deriving the correct and well-behaved master equa- 
tions in next-to-leading order perturbation theory. First 
we rewrite the equations by collecting the elements of 
the density operator into a vector, P, and the elements 
of the rate super-operators into matrices W, , L act- 
ing on this vector. Up to 4th order in the perturbation 



5 



expansion the equations can now be written as 

( -iL + W (2) + W (4) ) P, (19) 

(20) 



1 



e J P, 



I r = e T (w< 2) +W^)p. (21) 

The trace in Eq. $T$, JTHJ) is effected by the multipli- 
cation with the auxiliary vector e T = (1, . . . , 1, 0, . . . , 0) 
to sum up all vector elements corresponding to diago- 
nal density-matrix elements. The sum-rule on the kernel 
reads e T W = T . 



1. Elimination of Non-Diagonal Elements 

The crucial assumption for the following discussion is 
that the spectrum is free from accidental degeneracies in 
the following sense: all pairs of states a =/= b which have 
non-zero non-diagonal density matrix elements P a b are 
well separated in energy on the scale set by the tunnel- 
ing rates. Models for molecular transistors with discrete 
vibrational modes, such as the Anderson-Holstein model, 
satisfy this condition, provided that the vibrational level- 
spacing is larger than the tunneling coupling, since spin- 
selection rules generally prohibit coherence between the 
degenerate spin-states, unless broken by e.g. magnetic 
anisotropy or spin-polarization of the electrodes. One 
can always eliminate the non-diagonal elements and in- 
corporate their effect in a correction to the rates coupling 
diagonal elements. To this end we collect diagonal (d) 
and non-diagonal (n) density matrix elements into sep- 
arate vectors P. and P n , separate (|19[) into blocks and 
denote W = W( 2 ) + W( 4 »: 



~0d~ 




o„ 





P d 



(22) 



It is clear from (|15p that L is only non-zero in the nn 
block. The sum-rule implies 



a In 



(I ■ 

T 
n > 



(23) 
(24) 



where the multiplication with the vector ej = (1, . . . , 1) 
sums up all d- vector elements. We can now eliminate pro- 
cesses into the non-diagonal sector of the density-matrix 
by solving the equation from the lower block for the non- 
diagonal part of the density matrix, P n , and inserting 
this back into the equation in the upper block for the di- 
agonal part. Due to the clear separation of energy scales 
(non-degenerate spectrum) we can expand in the small 
quantity W„„L^. Consistently neglecting terms of or- 
der > H^, we then obtain an effective equation for P^: 



W d P d = 0, 



eJPd = 1, 



(25) 
(26) 



W rf 



W 



dd 



w 



(4) 
dd 



iW^L^wS- (27) 



The diagonal elements of the density matrix (vector of 
probabilities) thus satisfy what looks like a classical rate- 
equation, but with the effective rates (|27p . A completely 
analogous calculation for the correction to the current 
from non-diagonal elements gives 



£ (Wi, 



(W Jr ) 



w 



(2: 



(/</ 



w 



(4) 



dd 



(wg>) 



dn 



L" 1 W (2) 

%n vv rid ■ 



(28) 



(29) 



It can easily be shown that (|27|) and ([29]) are real, en- 
suring that the diagonal elements of the density matrix 
as well as the current are real. Due to Eq. (f2"51 12"4"|) the 
effective rate matrix satisfies the sumrule eJW d = T , 
so that Eq. (|25| with Eq. ([26]) determine the unique sta- 
tionary solution for the vector of diagonal density matrix 
elements (probabilities). Eq. (|25ti29[) form another cen- 
tral result of this work and we comment on their signifi- 
cance and importance. The advantage of the formulation 
in terms of effective rates, compared to solving Eq. (|19l 
[2~T]) directly, is threefold: (i) the effective rate matrices 
include the effects of coherences only up to order H^, 
just as the other effects of tunneling; (ii) it makes it ex- 
plicit that the 2nd order coherences effectively give 4th 
order effects in the rates for the occupations, something 
which is hidden in Eq. (TT9^) : (iii) it shows that the large 
matrix W„„, as well as all 4th order matrices which are 
not diagonal in initial and final state indices, need not be 
evaluated, significantly simplifying the calculation. The 
appearance of the correction in the effective rate has an 
intuitive physical meaning in the time-domain: it corre- 
sponds to a process starting (W nd ) and ending (W<j n ) in 
a diagonal state, through two tunnel processes. In the in- 
termediate non-diagonal state the free evolution involves 
rapid coherent oscillations at the Bohr-frequencies con- 
tained in L„„ (see (IT5T)). Due to the latter, these so-called 
non-secular terms [39j should be neglected in a lowest or- 
der approximation. However, these correction terms from 
coherences between non-degenerate states, although for- 
mally containing only 2nd order rates, contribute in 4th 
order to the occupancies, where they are crucial unless 
special model properties (conservation laws) make the 
matrix ~W dn vanish exactly. They scale in the same way 
as processes described by when one uniformly re- 

duces the tunneling matrix elements. Finally, we have 
found by numerical calculations for several model sys- 
tems that partial cancellations between the non-diagonal 
correction terms and diagonal 4th order terms are crucial 
for obtaining a physical result: if these corrections are ex- 
cluded one obtains SET-like resonances in the Coulomb- 
blockade regime below the inelastic cotunneling thresh- 
old. These are artifacts due to incorrect, large occupa- 
tions of the excited states, even at zero bias voltage. De- 
pending on the parameters of the model, negative occu- 
pation probabilities may even result, particularly when 
the tunneling amplitudes (JSJ7]) vary strongly from state 



to state. Accounting for the non-diagonal correction 
terms no such artifacts occur. Models for SMTs are typ- 
ical systems where the neglect of these non-diagonal cor- 
rections results in dramatic, spurious effects in the occu- 
pations and current. 

Summarizing: in the limit of large level-spacing, 
Eq. (|25If2T))) are the correct expressions for the occupa- 
tion probabilities and current. The corrections from 2nd 
order non-diagonal terms contribute only in 4th order 
in Ht- in a consistent 2nd order calculation they must 
be omitted whereas in a 4th order calculation they must 
be kept unless all non-diagonal elements vanish due to 
selection rules. 



2. Calculation of Diagonal Elements 

Having eliminated the non-diagonal elements, the re- 
maining problem is the solution of the kinetic equation 
for the diagonal elements ([25H26)) only. This requires 
some care since the effective rates (|27jl contain both 
2nd and 4th order terms, as was discussed in previous 
works {37l |40| (where corrections from non-diagonal el- 
ements were exactly zero due to selection rules). The 
problem is most easily understood from a simple exam- 
ple. Fig. [Ha) shows the result of 4th order perturbation 
theory for the single-level Anderson model in a magnetic 

field by solving Eq. (f25H29]) (due to spin-conservation 

(2) 

W^ n = and non-diagonal elements play no role) and 
Fig. [Hb) indicates relevant tunneling processes in re- 
gions (1), (2) and (3) in (a). In region (2) the excited 
spin-state can be populated by 4th order processes (in- 
elastic cotunneling). Since sequential tunneling out of 
this state is only possible at larger bias (region (3)), it 
can only relax by another inelastic cotunneling process 
back into the ground state. This latter process would 
not be included if one insists on an order-by-order solu- 
tion, i.e. expand also the occupation vector in powers of 
H T - P = <5P (0) + <5P (2) , solve for <*p(°) and 5P^ sepa- 
rately and discard the term W^SP^ in (j25|) . which is 
formally of order 6. Such an approach thus breaks down 
since the excited state is "pumped" up by inelastic tun- 
neling, but not allowed to relax by 4th order processes, 
yielding an unphysical solution: W^<5P(°) provides an 
inflow into the excited spin-state, but no outflow since 
jpC) i s only finite for the ground state, resulting in an 
artificially large correction SP^ 2 K Eq. ([23)) on the other 
hand has well behaved solutions, in which the occupancy 
of the excited spin-state is determined by the competi- 
tion of in- and out-going 4th order rates. We empha- 
size that the problem with the order-by-order solution is 
not of a numerical nature and occurs even if the equa- 
tions are solved analytically. It is of a general nature 
and occurs whenever all lowest order rates connected to 
some state are suppressed. Ref. (37j suggests dividing 
the Coulomb diamond into different regions, adapting 
the solution of the master equations thereafter, e.g. us- 
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FIG. 1: (Color online), (a): Occupation of excited spin-state, 
P+ (top, positive bias) and differential conductance (bottom, 
negative bias) for the single-level Anderson model where the 
spin degeneracy is lifted by an applied magnetic field. Here 
r L( r = r R(T = 10~ 2 T = 5 x 10~ 5 f7, where U is the charging 
energy, and e-f — ej = 50T. (b): Energy diagrams in the 
regions (1), (2) and (3), separated by green dashed lines in 
(a). In (1) only elastic cotunneling is energetically possible. 
In (2) also inelastic cotunneling can take place, but the only 
way for the excited system to return to the ground state is 
by another inelastic cotunneling process. In (3) the excited 
state can be emptied also by sequential tunneling processes 
(COSET). 

ing the order-by-order solution in the SET regime only. 
However, for a general SMT model such a division is not 
possible since even in the SET regime some rates may be 
suppressed by e.g. Franck-Condon or magnetic blockade 
effects. Always solving Eq. (|25p guarantees a physical 
solution in the sense that in- and out-going rates of all 
states are treated on an equal footing and the accuracy 
of the method is only limited by the order of the pertur- 
bation expansion of the kernel W. 

III. NON-EQUILIBRIUM 
ANDERSON-HOLSTEIN MODEL 

We now turn our focus to the Anderson-Holstein 
model, choosing the specific form of the molecular Hamil- 
tonian (JTJ 

H = e ]T dUa + jn(n - 1) + uj(tfb + i). (30) 

a 

The first two terms describe an electron in a single molec- 
ular orbital with electron operators dj., d a for spin a and 
h = d\d a denotes the number of excess electrons on 
the SMT. The last term describes the quantized vibration 
of the SMT through the operators b\ b. The cigenstates 
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FIG. 2: (Color online). Vibrational potentials in charge state 
N = and N = 1 and lowest corresponding vibrational wave- 
functions. The minimum of the potentials are shifted along 
the vibrational coordinate x by \/2\ (in units of the zero-point 
amplitude of the oscillator) . 

| a) in Eq. ((T|) thus have an electronic and a vibrational 
part, |o) = |e)|m e >, where |e> = |0>, | |>, I I), I U) 
denotes the electronic state with TV = 0, 1, 2 excess elec- 
trons on the molecule, and \m e ) labels the state of the 
oscillator. 

We have written Eq. (j30|) in the standard polaron-basis 
where e denotes the experimentally controllable effective 
energy level and U denotes the effective charging energy, 
both containing polaron-shift corrections [9(. The di- 
mensionless electron-vibration coupling, denoted by A, 
appears in an operator which displaces the vibrational 
states by V2X along the vibrational coordinate (normal- 
ized to the zero-point amplitude), whenever an electron 
tunnels from the electrodes onto the molecule, see Fig. [2] 
Thus the addition of an electron to the SMT induces a 
transition N — > N + 1, accompanied by a change of its 
vibrational state to' — > to. The matrix element for this 
process is reduced relative to the pure electronic tunnel- 
ing amplitude by the Franck-Condon overlap of vibra- 
tional wavefunctions in two different charge states of the 
SMT: 

/ w = (m| e - A ( fct - b V) (31) 

= ( - A )"' e -^^L™r™'(A 2 ), 

for to > ml (replace to <-> m' for to < to') where Lj(x) 
is the generalized Lagucrre polynomial. The dependence 
of the Franck-Condon factors on to and to' was discussed 
in detail in For a tunneling event starting from the 
vibrational ground state, to' = 0, the Franck-Condon fac- 
tors have a significant amplitude for A 2 — A < m < A 2 + A, 
corresponding to classically allowed transitions. More 
generally, for moderate to strong coupling there arc broad 
regions in the to, to' plane, bounded by the so-called 
Franck-Condon parabola, where vibration assisted tran- 
sitions have significant amplitude. These finite ampli- 
tudes for transitions to a range of vibrational states make 
the coherence between all pairs of these states important 
for the 4th order calculation, even though they are non- 
degenerate (on the scale of the tunneling broadening), 
i.e. there are many non-zero elements of Wd n and W nc ; 
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FIG. 3: (Color online). Differential conductance as a func- 
tion of gate and bias voltage close to the N — <-> N — 1 
degeneracy point, for A = 1 (a) and A = 2 (b). In the logarith- 
mic scale the lower end has to be chosen positive, preventing 
negative values from being correctly displayed, see e.g. white 
areas inside the sequential tunneling region in (a) which actu- 
ally correspond to very weak negative differential conductance 
(NDC). 

in (J22l). 

Here wc arc interested in transport close to a charge 
degeneracy point, accounting for the fact that the charg- 
ing energy together with the confinement-induced level- 
spacing typically constitute the largest energy scales in 
SMTs. We therefore restrict the model to electronic 
states with charge N = 0, 1, equivalent to taking U — ► oo 
in Eq. ([30]) . Without loss of generality we take e = — aVc, 
where a is the gate-coupling factor, i.e. we associate 
e = with zero gate voltage. The tunneling matrix el- 
ements for an electron tunneling onto the molecule are 
given by T^' + = S Szrjy /pt r f mm/ , where the eigenstates 
are labeled by the quantum numbers a = (s z ,m) in the 
N = 1 charge state and a' = (0, m') in the TV = 
charge state, with s z denoting the spin-projection of the 
molecule. We have everywhere used u> = 40T = 10 4 Fm, 
where Tm is the maximum sequential tunneling rate, i.e. 
r M = T x max(|/ mm ,| 2 ) and T = \2vJpt L \ 2 = \2iT^pt B \ 2 
is the pure electronic tunneling rate for symmetric cou- 
pling to the left and right electrodes. We set the width 
of the conduction band to D = 250uj. 



A. Intermediate Coupling 

The differential conductance as a function of gate and 
bias voltage is shown in Fig. [3] in the case of intermediate 
electron- vibration coupling, A = 1 in (a) and A = 2 in (b). 
The regimes where SET processes give the main contribu- 
tion to the current are triangle-shaped regions emanating 
from the point Vg = 0, where the energy for electron ad- 
dition without changing the vibrational quantum number 
lies between the electro-chemical potentials. Due to the 
quantized nature of the vibration of the SMT, additional 
sharp peaks appear in the differential conductance, asso- 



ciated with a change of the vibrational quantum number. 
Since in the SET-regime this is accompanied by a change 
in the charge, the positions of these peaks depend lin- 
early on the applied gate voltage. At the k-th resonance 
line (counting from V = 0) a new set of transitions be- 
comes energetically allowed, where the vibrational quan- 
tum number changes by k upon (dis)charging. 



Outside these two regimes, SET processes are sup- 
pressed by Coulomb blockade and one charge state is 
stable. Here no features are seen in a plot correspond- 
ing to Fig. [3] calculated to lowest non-vanishing order 
(not shown, see e.g. 0, SH). However, since we include 
all next-to-leading order processes, distinct features ap- 
pear in this region, which we now discuss. When the 
bias voltage reaches the vibrational level spacing, inelas- 
tic cotunneling processes exciting one vibrational quan- 
tum become energetically allowed. Due to the harmonic 
spectrum, this makes every excited vibrational state for 
fixed N accessible through a sequence of such tunneling 
processes: the molecular vibration is driven out of equi- 
librium. Each inelastic process involves the virtual oc- 
cupation of an adjacent charge state with an arbitrary 
vibrational excitation number. The onset of inelastic 
cotunneling is seen as steps in the differential conduc- 
tance, whose positions are independent of the gate volt- 
age since the process does not change the charge state of 
the SMT. The magnitude of the steps however depend 
on the gate voltage since the occupation of the virtual 
intermediate state is algebraically suppressed with the 
energy of this state. Similarly, at V = ku inelastic co- 
tunneling processes exciting k vibrational quanta become 
possible. The corresponding 2nd and 3rd inelastic co- 
tunneling steps are weakly seen for A = 2, while, for the 
tunneling coupling considered here, the suppression of 
the corresponding Franck-Condon factors renders them 
invisible for A = 1. 
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FIG. 4: (Color online), (a) and (b): Differential conductance 
as a function of gate and bias voltage for strong electron- 
vibration coupling, A = 3 in (a) and A = 4 in (b). For some 
values of the applied voltages the COSET processes result in 
closely spaced positive and negative differential conductance 
peaks, corresponding to peaks in the current, (c): Current as 
a function of bias voltage for A = 3 along the dashed green line 
in (a) , where the COSET processes give rise to a step + peak 
feature. As the vibrational relaxation rate, 7, is increased 
relative to Toi = r|/oi| 2 the peak vanishes, while the step 
remains. Inset: occupation of the vibrational ground state of 
the N = charge state, including also the result for A = 4 
without relaxation, (d): Sketch of lowest vibrational states in 
the N = 0, 1 charge states. An example of a COSET process 
contributing to the step + peak in (c) consists of an inelastic 
cotunneling process (blue arrows), followed by a sequential 
tunneling process (black arrow) into the vibrational ground 
state of the unstable charge state (N — 1). This may in turn 
sequentially relax (red arrow) to the vibrational ground state 
of the stable charge state (N — 0). 



A striking difference between Fig. [21a) and (b) is the 
appearance of gate-dependent lines inside the Coulomb 
blockade region in (b). The gate-dependence indicates 
that these lines are due to processes changing the charge 
state of the SMT, but they cannot be due to SET pro- 
cesses starting out from the vibrational ground state, 
since these are exponentially suppressed by Boltzmann 
factors (energy conservation). They originate instead 
from SET processes starting out from an excited vibra- 
tional state, which has previously been occupied by in- 
clastic cotunneling processes. This sequence of leading 
and next-to-leading order tunneling processes is called 
cotunneling-assisted sequential tunneling (COSET) [H, 
HE HE HE S3] j m the context of inelastic electron tunnel- 
ing spectroscopy (IETS) often referred to as phonon ab- 
sorption peaks, see Ref. [43[ and references therein. For 
even larger electron- vibration coupling these features be- 
come more pronounced as discussed in the next section. 



B. Cross-over to Strong Coupling 

The results of the calculations for larger electron- 
vibration coupling are shown in Fig. [H where A = 3 in 
(a) and A = 4 in (b). The most obvious consequence of 
a large electron-vibration coupling is the suppression of 
the low bias conductance. This Franck-Condon block- 
ade stems from exponentially vanishing overlap integrals 
(Franck-Condon factors) between low-lying vibrational 
states @,0,|45| which is seen in Fig. [Has a suppression of 
the degeneracy point peak (the differential conductance 
peak at V = Vq = 0). In the case of an equilibrium 
vibrational distribution and lowest order transport cal- 
culation, the current would increase exponentially with 
increasing bias voltage, until the blockade is lifted at 
around V/2 — (A 2 — X)u> = m corresponding to the first 
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large Franck-Condon factor f m o- However, when the vi- 
brational distribution is pushed out of equilibrium by se- 
quential tunneling processes, this significantly enhances 
the current compared to the case of equilibrium vibra- 
tions. Additionally, when next-to-leading order transport 
processes are taken into account, which was done using 
the Golden-Rule T-matrix approach in Ref. [|[ to study 
the strong coupling regime (A = 4 and A = 5), elastic 
and inelastic cotunncling processes change the exponen- 
tial suppression into an algebraic one. Cotunncling pro- 
cesses take place through high lying virtual intermediate 
vibrational states (m ~ A 2 ) which have a large overlap 
with the vibrational ground state, and the suppression 
of these processes is only algebraic with respect to the 
energy of the virtual state, and therefore with respect to 
A. 

The lowest inelastic cotunneling step is clearly seen for 
both A = 3 and A = 4. Additionally we find an anoma- 
lous signature of COSET processes. For low bias volt- 
age, just above the inelastic cotunncling threshold, these 
processes give rise to positive differential conductance 
(PDC) features, i.e. current steps, showing up as blue 
lines in Fig. 0J However, at larger bias voltages for A = 3 
we observe pairs of white and blue lines, corresponding 
to closely spaced lines of positive (PDC) and negative 
(NDC) differential conductance (see note on logscale in 
caption of Fig. [3]) . The nature of these line-pairs is more 
clearly seen in the current as a function of bias voltage, 
see red solid curve in Fig. 0Jc): COSET gives rise to a 
step in the current and, surprisingly, superimposed on 
it a peak. Such a peak has not been reported previ- 
ously to our knowledge and represents the central result 
of this section. We point out that all signatures in the 
current depend on a complicated interplay of a multitude 
of transport processes, also involving coherent superposi- 
tions of vibrational states. Basically, the peak arises due 
to a competition between leading and next-to-leading or- 
der transport processes and is closely related to the non- 
equilibrium vibrational distribution. This becomes clear 
from the inset of Fig. 0Jc) where we show the occupation 
of the vibrational ground state of the N = charge state 
for bias voltages around the peak. Although many vi- 
brational excitations are involved, the sketch in Fig.[4jd) 
gives an indication of the types of relevant tunneling pro- 
cesses. As the bias voltage exceeds the vibrational level 
spacing, inelastic cotunneling (blue arrows in Fig. H^d)) 
starts to deplete the ground state in favor of higher lay- 
ing vibrational states in the N = charge state (the first 
excited vibrational state acquires almost all of the prob- 
ability lost in the inset of (c)). Cotunncling processes 
starting from the excited states now give a significant 
contribution to the current which slowly increases with 
voltage. As one approaches the threshold for COSET 
from below, the current sharply increases as relaxation of 
these excited states into the N = 1 states by sequential 
tunneling (black arrow) becomes energetically allowed 
with large SET rates. If the FC-blockadc is not fully 
developed, a sequential tunneling process starting from 



N = 1 into the N = ground state may now follow with 
a larger rate than the inelastic cotunneling rate depicting 
the ground state, enhancing its occupation. As the volt- 
age moves through the COSET resonance this feedback 
increasingly suppresses the contributions from cotunncl- 
ing processes starting from excited states, thereby sup- 
pressing the current. As a result a thermally broadened 
peak occurs on top of the current step in Fig.QJc). More 
generally, such peaks appear when cotunncling processes 
start to become significant (A not too small) and compete 
with sequential tunneling processes, not fully suppressed 
by Franck-Condon blockade (A not too large). 

The effects of relaxation of the vibrational distribution 
due to a coupling to a dissipative environment, i.e. to 
substrate phonon modes, now has an interesting effect: 
as it is increased, at first it only suppresses the peak 
by disrupting the above competition. To illustrate this 
we have included a relaxation rate on a phcnomcnologi- 
cal level through an additional rate matrix W re i ax . This 
matrix is calculated in the same way as the tunneling 
rate matrix W, by performing an analogous perturba- 
tion expansion in the coupling to the dissipative bath, 7, 
with the difference that the bath operators are Bosonic 
rather than Fermionic. However, wc here restrict our- 
selves to the limit of weak coupling to the bath, 7 <C T, 
in which case we can stop this expansion at lowest non- 
vanishing order, analogous to Ref. and incorporate 
the result in the 4th order electronic rate matrix Wj^. 
We emphasize that such a simplified treatment becomes 
invalid as 7 ~ T since this requires treating coupling 
to the electron and phonon reservoirs on an equal foot- 
ing. The results for finite 7 is shown in the green dashed 
and blue dotted curves in Fig. 0Jc). The step only van- 
ishes when 7 > Toi = r|/oi| 2 (not shown), causing the 
first excited vibrational state to always relax before being 
emptied by sequential tunneling. The peak on the other 
hand depends on allowing several cotunncling processes 
to take place between relaxation events, and is thus much 
more sensitive to the coupling to the bath, thereby pro- 
viding an accurate experimental probe of the strength 
of the dissipative coupling. Additionally, since it only 
occurs within a range of A ~ 2 — 3 it also reveals infor- 
mation concerning the strength of the electron- vibration 
coupling. 

For A = 4, we find qualitatively similar results as pre- 
sented in Ref. Q . The degeneracy point is almost com- 
pletely invisible due to the strong Franck-Condon block- 
ade and the COSET processes do not give rise to peaks, 
but rather to PDC lines at low bias. The NDC lines seen 
at higher bias running perpendicular to the Coulomb di- 
amond edges occur already in a lowest order calculation 
within the sequential tunneling region 0. These NDC 
lines are seen to continue into the Coulomb blockade re- 
gion in our next-to-leading order calculation and are of 
a different origin. The absence of peaks at low bias is 
due to the fully developed Franck-Condon blockade, sup- 
pressing SET between vibrational ground states, thereby 
breaking the feedback mechanism which generates the 
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peaks. The pink, fine dotted line in the inset of Fig. [4|c) 
shows the ground state occupation for A = 4,7 = 0. It 
is clearly seen that, in contrast to the A = 3 case, the 
ground state does not become fully occupied above the 
threshold for COSET. 



IV. CONCLUSIONS AND OUTLOOK 

In this paper we have presented explicit kinetic equa- 
tions for quantum transport, valid for a generic class of 
molecular quantum-dot type systems, accounting for all 
contributions up to 4th order perturbation theory in the 
tunneling Hamiltonian and the complete non-equilibrium 
molecular density matrix. Due to the broadening of the 
states, which is treated correctly in the perturbation ex- 
pansion, all terms are automatically well-defined for any 
set of system parameters. The effective 4th order tran- 
sition rates, coupling diagonal elements of the molecular 
density matrix, include corrections from non-diagonal el- 
ements between non- degenerate states. In contrast to 
lowest order perturbation theory these corrections are 



essential for a physically correct solution. Applying the 
theory to the specific model of a molecular transistor cou- 
pled to a localized vibrational mode, we have shown that 
the signatures of cotunncling-assisted sequential tunnel- 
ing become more pronounced as the strength of the 
electron-vibration coupling is increased. In the cross-over 
to strong electron-vibration couplings, the cotunncling- 
assisted SET processes were shown to give rise to cur- 
rent peaks in the Coulomb blockade regime, which sig- 
nal a non-equilibrium vibrational state of the molecule. 
Their occurrence thus provides an indication of strength 
of the electron-vibration interaction. Since these peaks 
depend sensitively on an additional coupling to a dissi- 
pative bath, they also provide a way to experimentally 
estimate this coupling strength, 7, and thereby the im- 
portant Q-factor (Q = w/7). 
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from DFG SPP-1243, the NanoSci-ERA, the Helmholtz 
Foundation and the FZ-Jiilich (IFMIT). 



APPENDIX A: DERIVATION OF THE KINETIC EQUATION 



Our goal here is to derive the propagation of the reduced density matrix in Laplace space (jlip starting from 
Eq. (|10p . In the process we derive all diagrammatic rules. A number of techniques exist for calculating the trace 
over the reservoirs explicitly, such as projection operator techniques [46[ or path integral methods [47l. Although 
being formally equivalent to a diagrammatic expansion on a Keldysh double contour, see e.g. Ref. [33|], the diagram 
technique derived below has a number of advantages: (i) it is completely formulated and derived in Laplace space, (ii) 
a minimal number of diagrams represents all contributions in a given perturbation order, Keldysh and electron/hole 
indices being summed over, (iii) diagrams represent super-operators with diagram rules formally very similar to those 
for operators. This means we can postpone taking matrix elements, where the peculiarities of the Keldysh indices 
explicitly enter, to the end. Expanding the denominator in (|10[) we have 

P (z) = *Tr ( p + i -L T ^ -L T i + ....) P (0) p R , (Al) 

R [ z — Lr — L z — Lr — L z — Lr — L z — Lr — L J 

where (z — Lr — L)~ l is the free propagator and only even powers in Lx gi vc a non- vanishing contribution when 
performing the trace. The crucial step in developing a compact formalism is to ensure from the outset that Wick's 
theorem can be applied to super- operators in the same way as for operators. This is achieved by the definition of dot 
(G) and reservoir (J) super-operators by their action on an arbitrary operator A: 

a 2p £(N+pr)) 

rejT) P 2^ 2^ 1 r*( m ) \ -A|ai_)<a 2 -|, P=- ( ' 

N a lp <£N k 



= { °r^ A ' P Z + (A3) 

where we have assumed the tunneling matrix elements to be real- valued. Here p — ± is a Keldysh index, distinguishing 
between the forward (p = +) and backward (p = — ) time-evolution on a standard Keldysh double contour diagram. 
The index r\ = ± indicates an annihilation / a creation reservoir field operator. The product pr\ = ± has a physical 
meaning: when acting with G v rrjrt on a density operator, an electron / a hole is added to the dot from electrode r by 
projection between dot states with different charge and spin. The amplitude involves the tunneling matrix clement 
and a Keldysh sign p. An additional Keldysh sign p N ° appears in the amplitude. Importantly it can be assigned 
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FIG. 5: Diagrammatic representation of super-operator expressions. Processes evolve from right to left i.e. the diagrams have 
the same ordering as the expressions, (a): An example of an (irreducible) term in the expansion (|A6|I . (b): Separation into 
irreducible parts (self-energy or kernel, W (z)) and free evolution, II (z). The rightmost diagram is the only one contributing 
to the leading order self-energy, W^ 2 \ while the two other diagrams are the only ones in next-to-leading order, contributing to 



in any super-operator expression (i.e. without taking matrix elements) by simply counting the number Nq of Gs 
standing to the left (i.e. at later times). The explicit matrix elements of G (c.f. Eq. (TH]) ) required below are 



v ai + ai_ 



P 



1+N -N a 



(A4) 



where p = —p. Here the Keldysh sign is written as the parity of the charge difference between the final state of the 
G i.e. (— l) Ar ° 2 + ~ Na 2- — (_ i"j N ° (to see this, use that acting with Lt (~ G) changes the charge difference between 
the forward and backward contour of a Keldysh diagram by ±1, and that each diagram must start and end in a state 
which is diagonal in charge due to charge conservation of the total system). With these definitions it can be verified 
that the interaction Lt can be written as 



L T =^ 

prari 



dum Na G p J p 



Pi J i 



(A5) 



where in the second form we have defined the short-hand indices i = riair\iUJi and implicitly sum over Pi, ri,o~i, rji and 
integrate over Wj. The reservoir super-operators satisfy LrJ?* = Jf 1 (Lr — Xi) where Xi = rjitOi. In each term in the 
expansion, 



Tr- 



-Ln 



1 



R z — Lr — L z — Lr 



-Ln 



Pn 



.pf Gl [TtJP 
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■J? 1 PR 
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Z + Xr, 



L 



z — Lr — L 
1 



z + X n -i 



P(0)PR 
U n-1 ■ ■ ■ U 2 



1 G P1 —^—P (0) 

z + Xi - L 1 z-L V ' 



(A6) 



we can then pull all Js through to the left when adding Xi = x\ + X2 + . . . + Xi to Lr in the free propagators. 
Using Lr/»r = 0, pr can be pulled through as well. Since the reservoirs are assumed to be non-interacting we can 
now apply Wick's theorem to evaluate the trace over the super-operators J. In doing so one generates a Keldysh 



sign which exactly cancels p 



JV G „ 



N Gl 
■Pi 



This motivates including the canceling sign in the dot (|A2|) and tunneling 



Liouvillian (|A5|) super-operators to keep the final diagram rules simple. We contract pairs of reservoir super-operators, 
each contraction giving a factor 



la 



Pi(Jj 3 Lf )r = PiSr j r i 8 <T} <r i 6- fljtV .6(u)j - UJi)f(pi(Xi - 7?ijlt n )/T r< ), 



(A7) 



where f(x) = (e x + 1) _1 is the Fermi- function and T ri is the temperature of reservoir Ti (from hereon we assume equal 
temperatures of all reservoirs, T Ti = T). The Wick's sign follows in the usual way as the sign of the permutation 
which disentangles the contractions. All the Keldysh signs arise because the regular Wick's theorem can only be 
applied after all operators have been put on on the same forward Keldysh contour (i.e use cyclic invariance of the 
trace), see [3o| for details. Each super- operator in the expansion (|A6[) can thus be represented diagrammatically as 
usual by a directed free propagator line, [z + Xi — L) , interrupted by vertices G Vl which are contracted in pairs. 
A contraction of super-operators G^ 3 and Gf l with j > i is represented by an undirected line, see in Fig. [Ha). Since 
-rji,ujj = u>i,o-j = ai,rj = r 2 ; are enforced by the contraction, it can be unambiguously labeled by the indices 
(Tj, rji, r, of the earliest vertex G p ' . The sum in Xi collects only those x indices of lines passing over the free 
propagator segment i (contraction lines of one vertex to the left and one to the right), the other ones cancel. We now 
collect into W (z) all irreducible diagrams, i.e. those where any vertical cut will hit at least one contraction line, and 
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let II (z) = i(z — L) 1 be the contributions from free evolution of the molecule, see Fig. [5fb) . The molecular density 
matrix in Laplace space is now given by 



71 = V ' y ' 

where in the last step we have arrived at Eq. (fTTj) . The expectation value of the current operator I r is calculated 
analogously: 

(J r ) (z) = Trl r p (z) = TrL /r — P (0) p R = TrW Ir (z) P (z) . (A9) 

z — Lr — L — Lt m 

In contrast to Eq. (|10|) we trace over the full system, molecule + reservoirs (Tr = Trp(TrM)- Under the trace the action 
of the current operator I r on an arbitrary operator A has been expressed using the super-operator Lj r A = i |/ r , aX 
(aniz-commutator) which takes the same form as L^: 

Li r ^{G Ir )fJf\ (A10) 

Going through similar steps as above, we introduce a kernel Wj r (z) which differs from W (z) only by having the last 
G vertex replaced by a current vertex Gi r with matrix elements: 

((Gir) p r u m y i+ai ~ = {G^x: + a :; ■ ( An ) 

\ / 024-02- ^ 



APPENDIX B: DIAGRAMMATIC RULES AND PROPERTIES OF THE KERNEL 



The expression (| A8|) is still formally exact, but requires summing up all irreducible diagrams to obtain the kernel, 
which in general is not possible. We can write W (z) — W^ 2 ^ (z), where W^ 2k ^ (z) includes all terms with 

2k tunneling vertices, giving a pcrturbativc expansion in the tunneling Liouvillian Lt i.e. W^ 2k ^ ~ L T fc . We now 
summarize the diagrammatic rules obtained in Appendix [X] for calculating the zero- frequency z = iO contribution to 
the kernel: 

w*» w = e (jit) (-v N - G * i0+ xL-L G ™-i ■ ■■^icrk^i^- (B1) 

contr. 

Here one implicitly sums over all occurring Keldysh indices pi = ± as well as rj, <7j, r\i and integrates over all occurring 
energies Xi. 

1. (JIt) 1 Draw 2k vertices G 1 ^, i = 2k, . . . , 1 on a line. Connect pairs G v ^ 3 , G 1 ^ with j > i by a line denoting a 
Wick's-contraction. Equate the indices of G^ 3 to tw;, and — r)i and multiply by 

7 =Pif(Pi(%i - murJ/T). 

A vertex is contracted only to one other vertex and the contractions must be irreducible i.e. any vertical line 
through the diagram will cut at least one contraction line. 

2. (— 1) Determine the Wick's-contraction sign by counting the number of crossings of tunneling lines in the 
diagram. The parity of this number equals the parity of N p , the number of permutations required to disentangle 
the contractions. 

3. Assign a propagator (iO + Xi — L) _1 to segment i between vertex operators and G\ l . Here Xi = 
Si=coim xi is the sum of the energies of contractions passing through this segment i.e. the energies xi from 
all vertices i > I on the right contracted with some vertex to the left of the segment. 

4. X^contr • Perform 1-3 for every possible irreducible Wick's-contractions of the 2k vertices and sum them up. 

The current kernel Wj T is obtained by the same rules with the exception that the last vertex is replaced by the current 
vertex, G v ^£ — > {Gi^f-^k '• ^ ue ^° ^ ne additional 5-functions in the Gj r vertex (|A11[) . we need only include terms where 
an electron is added to the molecule from reservoir r in the final vertex (r/p = +). Additionally, due to the trace in 
Eq. (|A9[) we only need matrix elements which arc diagonal in final states. 
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One can check that these rules exactly reformulate the rules for the diagrammatic expansion of the kernels W and 
Wj r formulated previously (33], but in a compact manner well suited for constructing the general transport rates 
considered here. Fig. E^b) shows the single diagram for the leading order and the two diagrams making up 

the next-to-leading order kernel W^ 4 \ Since a Liouville diagram of order 2k has 2k Keldysh indices p, as well as k 
electron/hole indices rj, these diagrams account for 2^ x 2 rj = 8 (2 x 2^ x 2^ = 128) different Keldysh diagrams in 
leading (next-to-leading) order! The Keldysh representation is still useful in visualizing the character of the involved 
tunneling processes but need not be considered here. 

The general property of the kernel (W^z))^**"* - = (W(— z*))^~^ + * guarantees a Hermitian stationary-state 
density matrix. In the stationary limit [z — > iO), we have 

MW(iO))££: = Re(W(iO))£:£, (B2) 
hn(W(iO))2 + + 2Z = -MW(iO))^+. (B3) 

This has the important implication that elements of the kernel which arc diagonal in the double-indices ai+ = ai_ 
and a2+ = ai- are real-valued since they contain pairs of diagrams represented by complex conjugate expressions 
(obtained by inverting all p and r\ indices on a diagram). The same holds for Wj r (iO). 

Additionally, the charge difference between forward and backward Keldysh contours is conserved by each diagram. 
To see this, consider the action of a vertex operator G^ CTi7) . , which changes the charge number on contour pi by Pirji- 
Since in Eq. (|B 1[) this is contracted with G^ CT ._^. , this change in charge is either canceled (pj = pi) or equals that on 
the opposite contour (pj = —pi). The same hold for all pairs of contractions. 

APPENDIX C: 2ND ORDER 

In 2nd order, there is only one Liouville diagram, see Fig. [5^b). The diagrammatic rules give (we omit the argument 
iO) 

= -i l2l Gf— -Gf . (CI) 

tu + x\ — L 

The explicit evaluation of the matrix elements of this expression is discussed in some detail now, so that it can be 
skipped in the 4th order calculation where the expressions are less transparent, obscuring the basic simple operations. 
We introduce a shorthand notation for states on the forward / backward propagators: a, = a i+ aj_ and their energy 
difference E ai = E ai+ — E ai _ . Taking matrix elements and explicitly writing out summations and integrations, we 
obtain 

KOI - e e < g £->: (°s™): / ^ pl/( -o+'~'T )/r) < c2 > 

P2P1 ri(Ti?7i ai± ai 



V YYP2P1 X Vf 2B " lra ,T 01m ;° M s S a2 . (5, 

\ A^l ri<Tl(»7lp2) I"l<T 1 (77lpi) I "2p 2 "lp 2 



J -r 1 a 1 {rj 1 p 2 ) J 'r 1 <j 1 (,nPl) I " a 2p 2 a ip 2 " q 1pi a Oi»i 
P2P1 riiji ai± \ ai J 

x (-pi0((£? oi -m/iri)/r)-«v/(pi(£; oi -Wrj/T)), (C3) 

where 771 = —r\\ and p\ — —p\. The overall sign p\pi arises from several contributions. There is no Wick's sign since 
there is only one contraction (rule 2). The contraction-function gives a sign p\. Finally, the matrix elements of the 
vertices involve a sign pip\ and additionally a sign p\ since G^ 1 has an odd number of Gs standing to its left. 

For the integration we assume a flat density of states with a large bandwidth D 3> T,E ai — p r ,[i r — (i r i i.e. 
all energies E ai lie deep within this band, including all fj, r , meaning that we can neglect terms proportional to 

J D _ V dx^ «7/Dcl. Using x ^_ i0 = pi — inS (x), where P denotes the principal value, we split up the integral into 
real and imaginary parts. The imaginary part involves the Fermi-function and is the only contribution to elements 
of W diagonal in initial and final indices, which are just the well-known Golden Rule rates. The real part is only 
relevant for elements of W which are off-diagonal in initial or final indices, and involves the function (x rescaled by 
T) 

0(A) = -Re [ T dx- f ^ 



iO + x- X 



(] + 4) +ln 2^' ^ 
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where A = (E ai — i]ifi ri ) jT and ip is the digamma function. To arrive at this form we have used f(px) = (1 — p)/2 + 
pf(x) and neglected the integral Re ^|gq^rrf ^ XT/ D. Clearly, (f> (A) is symmetric for real- valued arguments, 

and we may write <p((E ai — rjifx rx ) /T) = (f>((i]iE ai — fx ri ) /T) i.e. only the distance of the addition energy to the 
Fermi-energy is relevant, irrespective of whether it is an electron / hole process (pi?7i = ±). The curve has a peak 
4>(0) = 7e + 21n2 + hijr^f = 1-96351 + ln^fy, where 7e is the Eulcr's constant, and logarithmic tails, <p(X) w ln^y for 
A> 1. 



APPENDIX D: 4TH ORDER 



In 4th order we have two irreducible contractions of the four vertices. We refer to the first diagram (leftmost in 
Fig. Efb)) as direct (D) type and the second one (middle in Fig. E](b)), which gets an additional sign from the Wick's 
contraction, as exchange (X) type. Applying the diagrammatic rules we obtain 

!*«<«) - -i^Cf aT i- T C5- - - - 1 - - ^ -^—^ (D) (Dl) 

+ ™" G " WT^-l g "' ,o + Xl l a - L G "l^7-L G '' ■ < x » ' D2 » 

Taking matrix elements, expanding all indices and explicitly writing out all summations and integrations this becomes 



P4P3P2P1 r 2 ri CT2<Ti r) 2 ni a, 3 ±a 2 ±a 1 ± 



< r Pi \«3 (r p 3 ( r p 2 \»I ( r pi \«o [[, , P2Plf(p2(x 2 - mVr 2 )/T)f(pi{xi - 7] lf l ri )/T) 

l^inaj a4 ^ 2l - 2CT2 j Q3 l^ 2 r 2CT J 02 l^xr 1( rj ai ^1^2 ( . Q + ^ _ ^ )( . Q + ^ + ^ _ + ^ _ ^ 

/^P4 \»8/^p3 \«2f™ 2 \«i C^Pl \ ao ff^^ P2Plf(p2(x 2 -r)2Vr 2 )/T)f( Pl (x 1 -T 1U l ri )/T) 
■ \yf) 2 r 2 a 2 ) ai ^iriffj 03 l°i) 2 r 2 <T 2 J Q2 lW n rifri J Ql / / ^1^2 



(iO + x 2 - E as )(i0 + xi + x 2 - E a2 )(i0 + xi — E ai 

(D3) 

Note that the two expressions differ only by the lower indices of vertex 3 and 4 and by the electron frequency x\ , x 2 
in the propagator connecting these vertices. 

For a non-degenerate spectrum, as discussed in the main text, we need only the expressions for with 04+ = 04- 
and ao+ = do-) which are guaranteed to be real- valued. We first give the final, explicit result, before discussing how 
to arrive there. 



Re 



[W^ j y ' y ' ^ ' ^ ' ^a4.p 4 a3p 4 ^a 3 p 3 a 2 p 3 ^a 2 p 2 aip 2 Saip 1 aop 1 P4Pl 



CI: I'Q ^ _ 

P4.P3P2P1 r 2 rt rt 2 r\\ a 3 ±a 2 ±a 1 ± 



P2P1 



Erpa 3p3 a 2p3 rpO, 2p2 ai P2 x ■> rpCi 4iPi a 3pi / j,ai P1 a 0pl \ 
r 2 o 2 (f) 2 p 3 ) r 2 <T 2 (r) 2 p 2 ) / j r\tj\{f\\pi) r\tj\(ji\p\) J 
o" 2 &i / 

F((E a2 - 77i/i ri - i]2Hr 2 ) /T, (E a3 - ?7i/i ri )/T) - F((£ a2 - ?/i^ ri - r) 2 p- r2 ) /T, (E ai - r)ifj, ri ) /T) 



+P1O--P2) 



{E aa ~E ai )/T 
F((E a3 - mf i ri ) /T) - F((£ 0l - TjrMn) /T)" 



(E a3 -E ai )/T 

\ rp a ip4, a 3 Pi rpO, 2p2 a lp2 X * rp a 3p 3 a 2 P3 rp a lpi a O P1 \ 

' ' r 2 cr 2 (fj 2 p 3 ) r 2 a 2 (-q 2 p 2 ) / , riai(fjip4.) r 1 a 1 (r] 1 p 1 ) I £* 2 Pl 



\ C 2 (Ti / 

F (( E a 2 ~ »?lMri ~ mVr 2 ) /T, (E ai - ?/lMn) /T) - F((£? aa + Ki - JJiMn - ??2Mr 2 ) - J]l fi ri ) /T) 



(E a2 - E a3 - E ai ) /T 

F((E a2 - - ??2/ir 2 ) /T, (gag - ??2Mr 2 ) /T) - F((E a3 + £», - TJlpiri - T\2^r 2 ) /T, (Eg 3 - 7? 2 /X r2 ) jT) 

(Ea 2 — E a3 — E ai ) jT 



(D4) 
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where only two types of functions enter 

F(\',\) = 7r{0(A'-A)/(A) + &(A')[0(A'-A)-0(-A)]} (D5) 

(-A)/(A)-^0(-A)| for A'^0, 
F(X) = |0(A), (D6) 

where /(A) = (e A + 1) _1 and 6(A) = (e A — 1) _1 are the Fermi- and Bosc-function respectively and 4>(X) is given by 
Eq. (|C4[) . All expressions arising from the integrals are explicitly seen to be well behaved, since they take the form of 
differential quotients: whenever a denominator vanishes, the numerator also vanishes with the same power, resulting 
in a finite value. The rates are thus well-behaved functions of all model parameters including the voltages. 

We now discuss the steps leading from Eq. (|D3[) to Eq. (|D4[) . The tunnel matrix elements enter automatically via 
the vertices Eq. (|A4[) . The four vertices give a sign P4P3P2P1, and the vertices G3 3 and G^ 1 give an additional sign 
PsPi (since they are followed by an odd number of vertices towards the left). Combined with the contraction signs 
P2P1 we get in total a sign p^pi for both diagrams. 

The remaining task is to obtain the closed- form expressions for the imaginary part of the two integrals. Normalizing 
the integration variables to T and then shifting them introduces the energy denominators Ai = (E ai — rjip ri ) /T 
and A2 = (E a2 — r\\ii ri — ?72Mr 2 ) /T. For the last propagator we get A3 = (E a3 — ?7i/i ri ) /T for the D type and 
A3 = (E a3 — T]2fJ-r 2 ) /T for the X type diagram. The integrals are then split into partial-fractions: 



D TJJ A 3 -Ax \iO + x 1 + x 2 -\ 2 ) yiO + X! -A 3 iO + n-Ai 



iO + xi — Ai iO + X2 — A3 



TJJ A 2 — A 3 — Ai \i0 + x\ + x 2 — A2 iO + x\ + x 2 - A 3 - Ai 

1 ^ (D8) 

where Iq denotes the integral in (|D3j) in the D type and Ix the one in the X type diagram. These can be expressed 
in the integrals encountered in 2nd order. This is done most efficiently by first noting a number of sumrules which 
are satisfied by the integrals (but not by the diagrams!) in the wide-band limit: 

£ I£T = £ ir 1 - E lP x Pl = °- (D9) 

Pl=± Pl=± P2=± 

Summing the integrals over a Kcldysh index we eliminate one Fermi- function using Y] ± / (piXi) = 1. We can 

then first evaluate the integral over a;^ on the same contour as for the 2nd order integral, see Appendix |EJ If the 
integrand vanishes faster than x~ x the contribution can be neglected in the wide-band limit, even when performing 
also the 2nd integral. From the original expressions for the integrals in Eq. (|D3|) one sees that this is the case, except 
for the integrand Id considered as function of x 2 - Therefore J2 P2 =±^D Pl ^ This implies that Ix is proportional 
to P1P2 , while Id additionally contains a term proportional only to p\ : 

_ 1 F(X^)-F M) + 1„ A)-^ (m()) 



T A3 — Ai T A3 — Ai 

rP2Pl 1 F(A 2 , AQ - F (A 3 + Ai, Ai) + F(A 2 , A 3 ) - F(A 3 + A 1; A 3 ) 

Ix ~ T P2Pl A 2 -A 3 -A! ' (DU) 



It remains to be shown that ^(A', A) and F (A) actually are given by Eq. ()D5|> and ()D6P respectively. To do this we 
now use the expansion 

/ (P2X2) f (pixi) = P2P\! {xi) f (x 2 ) + -pi (1 - pa) / + \p2 (1 - Pi) f [x 2 ) + const. (D12) 

As was noted when deriving the above sumrule, only terms containing the product / [x\) f (x 2 ) give a non- vanishing 
contribution to the X-type integral, and we only have to consider integrals of the form 

f{x')f{x) 



F(A',A) = Im [[ dxdx'- 

v ' J J (iO + x + x' -\')(iO + x-X) 



— 7rRc 



iO + x' - A' + A w J iO + x - A 

tt {<f>(\' - A)/(A) + b(X') [<f>(\' - A) - cj>(-\)}} . (D13) 
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Here we expanded Im {{x + iO) (y + iO)) 1 = — 7rRc |<5 (x) (y + iO) 1 + S (y) (x + iO) x | and used the relation 

/ (a/) / (x) = (/ (— x') — f (x)) b(x + x') in the second term. 

The D-type integral gives a non- vanishing contribution also for terms containing only f (x\). This yields the 
additional integral where / (P2X2) / (pi^i) \v\ (1 — P2) f 

F (A) = hm [[ dxdx'- 



2 J J (iO + x + x' — A') (iO + x — A) 

= -7r-Re / dx - = -<j> (A) . (D14) 

2 y iO + x-A 2 VK ' y ' 

Note that A' drops out of the answer since Rc J D ^ T dx' i0+x + x ,_ x , van i sncs f° r D/T ^> A, A'. 



APPENDIX E: CONTRACTION INTEGRAL 



We comment on the calculation of the integral: 



/ 



dx- 



iO + x - A 



= +Re-0 



1 A 

2 +Z 27 



-In- 



2^T 



i7T/(A) 



(El) 



It can be calculated with a smooth Lorentzian cutoff of width D/T and the result must then be expanded in the 
small parameter XT/D to lowest order [48| . This however involves unnecessary complications since the energy scale 
separation D/T ^> A is only used at the end. Here we indicate how this may be avoided, simplifying this and other 
similar calculations. We first note that although — iirf (A) clearly stems from ilm z _^ +iQ = —iir6 (z — A) one should 
not separate real and imaginary parts until the end of the calculation. We apply the residue theorem for a contour 
along the real axis and finite semi-circle in the upper half-plane i.e. not containing the pole z = A — iO. We obtain 
the integral with a sharp cutoff: 



D/T 

dx- 

-D/T X — \ + l\) 



dipz 



/(*) 



z - A + iO 



z=De i 'P/T 



fe=0 



z - A + iO 



(E2) 



z=iir(l+2k) 



where ko = [^pf — l\ (M denotes the integer part). We now explicitly calculate the contribution to the contour 
accounting for D/T ^> A. The latter is trivial since / (z) is equal to 1 for ? < argz < ir and elsewhere for z on 
a semi-circle of radius D/T 3> A, as one easily verifies. Since the remaining part of the integrand is independent of 
argz on this contour, we get a contribution —if - In the limit D/T ^> A the summation over Matsubara-poles can be 
extended to infinity and gives a digamma function plus a term depending logarithmically on the band-width: 




lnfc 



D 



1 1 .AN . D 

-ib - +i — +ln 

1 2 2tt 2-kT 



2' 



(E3) 



where we added and subtracted the Euler's constant, 7e = linin^oo Y^k=i 1/fc ~ Inn. The contribution from the arc 
cancels part of the imaginary part Inr0(l/2 + ix) = tt tanh(7ra;) /2 = 7r(l/2 — f(2irx)). In contrast, if one takes a cutoff 
function to make the integral vanish along the semi-circle for infinite radius (48[, one unnecessarily complicates the 
evaluation of the residues. 
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